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Abstract 

A  detailed  steady- state  isothermal  two-dimensional  model  of  a  proton  exchange  membrane  fuel  cell  has  been  developed.  A  finite  element 
method  was  used  to  solve  this  multi-component  transport  model  coupled  with  flow  in  porous  medium,  charge  balance,  electrochemical  kinetics, 
and  a  rigorous  water  balance  in  the  membrane.  The  model-predicted  fuel  cell  performance  curves  are  compared  with  published  experimental 
results  and  a  good  agreement  was  found.  The  complex  water  balance  in  the  membrane  was  investigated  and  the  operating  conditions  where  the 
membrane  becomes  dehydrated  were  identified.  The  effects  of  channel  width  and  bipolar  plate  shoulder  dimensions,  porosity,  and  the  relative 
humidity  of  the  inlet  streams  on  the  fuel  cell  performance  are  evaluated.  It  was  found  that  smaller  width  channels  and  bipolar  plate  shoulders 
were  required  for  high  current  density  operations.  As  the  electrode  area  under  the  bipolar  plate  shoulder  increases,  the  fuel  cell  benefits  more 
from  higher  porosity  electrodes.  The  anode  gas  stream’s  relative  humidity  was  found  to  be  more  critical  for  fuel  cell  performance  than  the 
cathode  gas  relative  humidity. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Since  the  early  1990s  there  has  been  growing  interest 
in  modeling  proton  exchange  membrane  (PEM)  fuel  cells. 
Earlier  models  by  Bernardi  and  Vebrunge  [1]  and  Springer 
et  al.  [2]  are  fundamental  studies  toward  understanding 
PEM  fuel  cells.  The  two  groups  both  developed  a  one¬ 
dimensional,  isothermal  model  of  a  membrane  electrode  as¬ 
sembly  (ME A).  Bernardi  and  Vebrunge  [1]  studied  the  net 
water  flux  through  the  membrane  and  Springer  et  al.  [2]  in¬ 
cluded  variable  membrane  hydration  in  his  model.  Fuller  and 
Newman  [3]  developed  a  two-dimensional  (2D)  model  to 
study  water  and  thermal  management  issues  and  used  con¬ 
centrated  solutions  theory  for  transport  in  the  membrane. 
Nguyen  and  White  [4]  studied  the  effects  of  various  forms 
of  gas  humidification  on  cell  performance,  heat  and  water 
management  and  reported  that  back  diffusion  of  water  from 
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the  cathode  to  anode  is  insufficient  to  keep  the  membrane 
hydrated  at  high  current  density  operations.  They  conse¬ 
quently  concluded  that  the  fuel  and  oxidant  must  both  be 
humidified. 

Variable  degrees  of  water  flooding  in  the  catalyst  layers 
and  electrode  backing  regions  was  studied  by  Baschuk  and 
Li  [5]  in  a  one-dimensional  model.  Later  Baschuk  and  Li  [6] 
also  examined  the  issue  of  carbon  monoxide  poisoning  in 
PEM  fuel  cells. 

In  recent  years,  a  general  trend  of  using  computational 
fluid  dynamics  (CFD)  to  model  PEM  fuel  cells  has  evolved. 
Gurau  et  al.  [7]  developed  the  first  real  two-dimensional 
model  of  a  fuel  cell  with  flow  channels  and  ME  A.  This 
“along- the-channel”  model  studied  the  effects  of  composi¬ 
tion  change  of  the  reactants  inside  the  channels.  Um  et  al.  [8] 
developed  a  two-dimensional  transient,  “along- the-channel” 
model  and  studied  the  change  of  current  density  with  chang¬ 
ing  cell  potential.  Um  and  Wang  [9]  extended  the  work  to  the 
third  dimension  and  also  studied  the  effects  of  flow  channel 
geometry  and  layout.  A  group  at  the  Electrochemical  Re- 
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Nomenclature 

a 

water  activity 

c 

concentration  (mol  m-3) 

cw 

mass  concentration  of  water  in  the  membrane 
(kgm-3) 

Dij 

binary  diffusivity  (m2  s-1) 

£>w 

water  diffusivity  in  the  membrane  (m2  s-1) 

EWm 

equivalent  molecular  weight 

F 

Faraday’s  constant  (96487  C  (mol)-1) 

I 

local  current  density  vector  (A  m-2) 

k 

exchange  current  density 

k 

constant  reaction  parameter 

kp 

permeability  (m2) 

L 

length  (m) 

nd 

drag  coefficient 

M 

molecular  weight  (kg  mol-1) 

N 

molar  flux  vector  (mol  m-2  s-1) 

Vv 

net  water  mass  flux  vector  (kgm-2  s-1) 

P 

pressure  (Pa) 

R 

universal  gas  constant  (8.314  Jmol-1  K-1) 

T 

temperature  (K) 

u 

velocity  vector  (ms-1) 

Eoc 

open  cell  potential  (V) 

w 

weight  fraction 

X 

mole  fraction 

Greek  letters 

a 

transfer  coefficient 

a 

conductivity  (S  m-1) 

£ 

porosity 

y 

reaction  order 

X 

membrane  water  content  (mol  H2O  (mol 
SO3-)-1) 

p 

viscosity  (kgm-1  s-1) 

1 1 

overpotential  (V) 

p 

density  (kgm-3) 

Subscripts 

avg 

average  value 

Uj 

components,  H2,  and  H2O  for  the  anode,  O2, 
H2O,  and  N2  for  the  cathode 

s 

gas  distribution  electrode 

m 

membrane 

w 

water  in  the  membrane 

Superscripts 

eff 

effective 

ref 

reference  conditions 

search  Center  of  Pennsylvania  State  University  developed  a 
large-scale  CFD  model  [10]  and  studied  the  two-phase  trans¬ 
port  issues  in  PEM  fuel  cells  [11].  Berning  et  al.  [12]  devel¬ 
oped  a  three-dimensional  model,  conducted  parametric  stud¬ 
ies  on  operating  pressure  and  temperature  as  well  as  geomet¬ 


rical  and  material  properties  [13],  and  developed  a  two-phase 
model  [14]. 

All  of  the  CFD  work  done  by  the  previously  cited  authors 
was  aimed  toward  understanding  the  PEM  fuel  cell  and  sim¬ 
ulating  the  performance  of  a  given  design.  However,  those 
models  were  not  suitable  for  optimization  of  the  fuel  cell  de¬ 
sign  and  other  units  that  comprise  the  fuel  cell  system.  To  our 
knowledge  optimization  of  PEM  fuel  cell  design  using  CFD 
has  been  attempted  only  by  Grujicic  and  Chittajallu  [15]. 
However,  the  model  used  by  Grujicic  and  Chittajallu  [15] 
assumes  constant  membrane  hydration  meaning  that  mem¬ 
brane  conductivity  changes  and  the  complex  water  balance 
was  not  investigated. 

Currently,  there  are  several  commercial  CFD  packages 
from  various  vendors  that  can  be  used  for  modeling  PEM 
fuel  cells  [16-19].  The  nature  of  modeling  PEM  fuel  cells 
requires  uncommon  water  balance  equations,  which  these 
commercial  CFD  packages  can  add  and  solve.  CFD  mod¬ 
eling  is  computationally  demanding  especially  for  three- 
dimensional  models.  Because  the  geometric  dimensions  of 
the  various  fuel  cell  regions  range  several  orders  of  mag¬ 
nitude,  the  models  require  very  small  elements  to  capture 
the  details  with  accurate  resolution.  For  example,  the  typi¬ 
cal  anode  and  cathode  catalyst  layer  thickness  of  PEM  fuel 
cells  are  on  the  order  of  10-20  jam,  however,  the  height  and 
width  are  10-20  cm  requiring  at  least  a  million  elements  for 
even  a  small  7  cm  by  1  cm  section  of  a  fuel  cell  [10].  The 
computational  cost  of  CFD  modeling  an  entire  cell  is  large 
and  we  are  not  yet  able  to  simulate  an  entire  fuel  cell  stack. 
However,  the  interactions  between  cells  in  a  stack  are  critical 
for  water  and  thermal  management  since  the  cells  located  at 
the  middle  section  of  the  stack  will  require  extra  means  of 
cooling. 

Considering  the  computational  cost  of  CFD  modeling  of 
PEM  fuel  cells,  using  a  three-dimensional  fuel  cell  stack 
model  is  not  practical  in  most  computing  environments.  The 
purpose  of  this  study  is  to  develop  a  fast,  robust,  and  de¬ 
tailed  two-dimensional  CFD  model  that  can  predict  effects 
of  channel  geometry,  and  water  management,  and  that  can 
be  used  for  optimization  of  the  fuel  cell  components,  design, 
and  operating  conditions  for  different  applications. 


2.  Problem  description 

A  schematic  illustration  of  a  PEM  fuel  cell  divided  into 
seven  sub-regions  is  shown  in  Fig.  1.  The  membrane  elec¬ 
trode  assembly  consists  of  an  ion  exchange  membrane  with  a 
thickness  of  100-250  p,m  sandwiched  between  two  Pt  based 
catalyst  layers  with  a  thickness  of  10-20  |jim  and  backed  by 
porous  gas  distribution  electrodes  (GDE),  which  is  usually 
carbon  cloth  with  a  thickness  of  200-300  pan.  The  MEA  is 
then  sandwiched  between  two  bipolar  plates.  These  plates 
serve  as  a  manifold  for  the  transfer  of  feed  and  products  in 
and  out  of  the  fuel  cell  and  as  the  current  collector. 


G.H.  Guvelioglu,  H.G.  Stenger  /  Journal  of  Power  Sources  147  (2005)  95-106 


97 


e 

» 

LOAD 


Bipolar  Plate 


Gas-Diffusion  Electrode 

Catalyst  Layer 

Membrane 


Bipolar  Plate 


Gas-Diffusion  Electrode 
Catalyst  Layer 


Anode 


Cathode 


Fig.  1.  Schematic  illustration  of  the  PEM  fuel  cell. 


The  heart  of  a  fuel  cell  is  the  active  catalyst  layers,  since  at 
these  regions  the  four  ingredients:  membrane,  catalyst,  elec¬ 
tronic  conductor,  and  chemical  reactants,  required  for  the  fuel 
cell  electrochemistry,  exist  together. 

Hydrogen,  which  is  usually  humidified,  enters  the  anode 
gas  chamber,  transports  through  the  porous  electrode  by  con¬ 
vection  and  diffusion,  and  dissolves  into  the  membrane  phase 
of  the  anode  catalyst  layer.  The  hydrogen  is  subsequently  ox¬ 
idized,  the  membrane  is  replenished  with  protons,  and  the 
carbon  conductor  receives  electrons  by  the  anode  reaction: 

2H2  — >  4H+  +  4e“  (R-l) 

The  oxidant,  O2  in  air,  which  is  usually  humidified,  enters  the 
cathode  gas  chamber,  transports  through  the  porous  electrode 
by  convection  and  diffusion  and  dissolves  into  the  membrane 
phase  of  the  active  catalyst  layer.  Protons  in  the  membrane 
pores,  coming  from  the  anode  side,  react  with  dissolved  O2 
at  catalyst  sites  in  the  active  layer  to  produce  water  in  the 
electrochemical  reaction: 

02  +  4H+  +  4e“  -*  2H20  (R-2) 

A  comprehensive  review  of  PEM  fuel  cell  design  and  manu¬ 
facturing  is  done  by  Mehta  and  Cooper  [20]. 

The  design  of  PEM  fuel  cells  requires  a  strong  understand¬ 
ing  of  the  processes  such  as  mass,  momentum  and  energy 
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Fig.  2.  Computational  domain  and  governing  equations. 
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transport,  electrochemical  reactions,  and  charge  balance  tak¬ 
ing  place  inside  the  fuel  cell.  A  PEM  fuel  cell’s  performance 
is  significantly  affected  by  the  water  content  of  the  mem¬ 
brane  due  to  the  nature  of  the  membrane  material.  Greater 
membrane  hydration  means  higher  conductivity,  less  poten¬ 
tial  loss  at  the  membrane  thus  higher  cell  voltage  and  power. 
Water  balance  in  the  fuel  cell  is  a  complex  phenomenon  and 
needs  to  be  controlled  accurately  for  the  successful  design 
and  operation  of  fuel  cells. 

A  typical  cell  thickness  is  in  the  range  of  1-2  mm,  this 
includes  the  seven  sub-regions  as  well  as  the  bipolar  plates. 
The  length  and  width  of  a  typical  single  cell  is  between  10 
and  20  cm.  This  order  of  magnitude  difference  in  dimensions 
requires  simplification  of  the  fuel  cell  to  2D  geometry  in  order 
to  minimize  the  number  of  mesh  elements.  The  thickness  of 
the  active  catalyst  layers  is  about  10  pan,  which  is  very  small 
with  respect  to  other  cell  components,  and  can  be  treated  as 
reactive  boundaries.  The  selected  2D  domain  includes  the  gas 
distribution  electrodes,  catalyst  layers  as  reactive  boundaries 
and  membrane  and  shown  in  Fig.  2. 

Model  assumptions: 

•  steady- state  operation, 

•  isothermal  operation, 

•  ideal  gas  mixtures, 

•  single-phase  model, 

•  isotropic  and  homogeneous  electrodes  and  membrane, 

•  the  membrane  is  considered  impermeable  for  the  gas 
phase  which  was  supported  by  the  findings  of  an  earlier 
study  [1], 

•  negligible  contact  resistance, 

•  minimal  membrane  swelling. 

Model  simplifications: 

•  two-dimensional  model, 

•  catalyst  layers  as  reactive  boundaries. 

The  PEM  fuel  cell  model  developed  is  a  comprehensive 
two-dimensional,  isothermal,  steady-state  model  providing  a 
detailed  description  of  the  following  transport  phenomena: 


3.  Model  equations 

3.1.  Gas  distribution  electrodes 

In  the  gas  distribution  electrodes,  Darcy’s  law  is  used  to 
model  the  flow  in  this  porous  media  with  the  pressure  gradient 
as  the  driving  force.  In  a  porous  structure,  the  global  trans¬ 
port  of  momentum  by  shear  stresses  in  the  fluid  is  negligible 
because  the  pore  walls  impede  transport  of  this  momentum 
to  the  fluid  outside  the  individual  pores.  Since  a  detailed  de¬ 
scription  at  the  resolution  of  a  pore  is  not  practical  in  most 
models,  homogenization  of  the  porous  and  fluid  media  is  a 
common  approach.  Darcy’s  law  is  based  upon  homogeniza¬ 
tion  of  the  porous  and  fluid  media  into  one  single  medium 
and  does  not  require  a  detailed  geometrical  description  of  the 
pore  structure. 

Darcy’s  law  states  that  the  velocity  vector  is  determined 
by  the  pressure  gradient,  the  fluid  viscosity,  and  the  structure 
of  the  porous  media  represented  with  the  following  equation: 

kr\ 

u  = - v-Vp  (1) 

where  u  is  the  velocity  vector,  kp  the  permeability,  /x  the  gas 
viscosity,  and  p  is  the  pressure.  Variables,  descriptions,  and 
their  units  are  given  in  Nomenclature  section  of  this  paper. 

For  multi-component  diffusion  in  gases  at  low  density  it 
has  been  shown  that  Maxwell-Stefan  is  a  good  approximation 
[21]. 

N  i 

V*;  =  -  V]  ——(xjNj  -  xtNj ),  i  =  1,  2, . . . ,  N  (2) 

pt  cDU 

where  D(j  is  the  binary  diffusivity  of  i  and  j,  c  the  concentra¬ 
tion,  Xi  the  mole  fraction  of  component  /,  and  N[  is  the  molar 
flux  vector  of  component  i.  Experimentally  obtained  binary 
diffusivities,  D?-  at  atmospheric  pressure  j7a tm  and  reference 
temperature  Tq  shown  in  Table  1  are  scaled  to  operating  tem¬ 
perature  and  pressure  according  to  [21]: 

D'l  =  ^AT«- P«>J  (jr)  ® 


•  multi-component  flow, 

•  diffusion  of  reactants  through  the  porous  electrodes, 

•  electrochemical  reactions, 

•  transport  of  electrons  through  the  electrodes, 

•  water  balance  in  the  membrane. 


The  species  balances  in  the  porous  gas  diffusion  electrodes 
are  solved  with  the  following  equation: 


(4) 


The  equations  governing  these  processes  include: 

•  ionic  balance  in  electrodes  and  membrane, 

•  the  Maxwell-Stefan  equations  for  multi-component  dif¬ 
fusion  and  convection  in  gas  distribution  channels  and 
gas  distribution  layers, 

•  Darcy’s  law  for  the  flow  of  species  in  porous  electrodes, 

•  water  balance  and  water  flux  in  the  membrane  governed 
by  diffusion,  convection,  and  electro-osmotic  drag. 


Table  1 


Binary  diffusivities  and  reference  temperatures  at  1  atm  [13] 


Gas  pair 

Reference 

Binary  diffusivity, 

temperature,  To  [K] 

Dij(T0,po )  [m2  s-1] 

D° 

^h2-h2o 

307.1 

9.15  x  1CT5 

D° 

^o2-h2o 

308.1 

2.82  x  1CT5 

D° 

^o2-n2 

293.2 

2.2  x  10-5 

D° 

^h2o-n2 

307.5 

2.56  x  IQ"5 
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where  i  =  H2  and  H2O  for  the  anode  and  i  =  O2,  H2O,  and  N2 
for  the  cathode  side,  w  the  weight  fraction,  and  p  is  the  gas 
mixture  density  calculated  with: 


£/*/  •  MW/ 


RT 


However,  due  to  the  porous  structure  of  the  electrodes  the 
binary  diffusivities  need  to  be  corrected  for  the  porosity  s 
of  the  electrode  media.  This  is  done  with  the  Bruggeman 
correlation  [1]: 


Dus 


1.5 


where  Djj  are  the  binary  diffusivities  from  Table  1,  corrected 
for  the  operating  temperature  and  pressure  with  Eq.  (3). 

The  Maxwell-Stefan  and  Darcy’s  equations  are  coupled 
through  the  velocity  vector  m,  and  density  p.  In  addition  to 
mass  and  momentum  balances  in  the  porous  electrodes  the 
charge  balance  is  also  solved  in  the  electrodes  with  the  con¬ 
tinuity  relationship: 

-V  •  (aseffV0s)  =  0  (7) 

where  of ff  is  the  effective  conductivity  and  </>s  is  the  potential 
in  the  electrode. 


3.2.  Membrane 

The  permeability  of  the  membrane  to  hydrogen,  oxygen, 
and  nitrogen  is  low  and  can  be  neglected  [  1] .  Thus,  only  water 
and  protons  are  transported  and  each  obeys  the  principle  of 
mass  conservation,  with  the  charge  balance  in  the  membrane 
accounted  for  with  the  continuity  relationship: 


-V  •  (ormV0m)  =  0 


where  0m  is  the  membrane  potential  and  crm  is  the  mem¬ 
brane  conductivity  which  is  a  function  of  membrane  water 
content  A,  i.e.,  the  local  [H20]/[S03_]  ratio  in  the  mem¬ 
brane.  Springer  et  al.  [2]  presented  the  following  empirical 
expression  for  the  conductivity  of  Nafion®  membrane  based 
on  their  experiments: 


or30  =  0.5139A  —  0.326,  for  A  >  1 


where  <730  is  the  membrane  conductivity  [Sm  ]]  measured 
at  30  °C  and  then  corrected  to  operating  temperature  T  with: 


<rm  =  exp 


(10) 


At  steady-state,  the  water  balance  in  the  membrane  re¬ 
duces  to: 


V  •  Nw  =  0 


(11) 


The  mass  flux  of  water  in  the  membrane  is  governed  by 
three  processes:  electro-osmotic  drag  due  to  the  flow  of  pro¬ 
tons,  diffusion  due  to  a  concentration  gradient,  and  convec¬ 
tion  due  to  pressure  gradient.  These  three  processes  lead  to: 

Hd  •  /  •  Mh20  _  „ 

Nw  —  Dw  V  cw  cwuw  (12) 

F 

where  cw  is  the  mass  concentration  of  water,  Dw  the  diffu¬ 
sion  coefficient  of  water  in  the  membrane,  nd  the  electro- 
osmotic  drag  coefficient,  I  the  local  current  density  vector, 
Mh2o  the  molecular  weight  of  water,  and  F  is  the  Faraday’s 
constant.  The  velocity  of  the  water  in  the  membrane,  uw  fol¬ 
lows  Darcy’s  law  Eq.  (1). 

The  membrane  water  diffusivity  is  related  to  the  temper¬ 
ature  and  water  content  of  the  membrane  with  the  following 
relationship  of  Motupally  et  al.  [22]. 


Dw 


3.1  x  10  7A[exp(0.28A)  —  1]  exp(— 2346/T),  A  <  3 

4.17  x  1(T8A[1  +  161  exp(— A)]  exp(— 2346/  T),  A  >  3 

(13) 


The  drag  coefficient  of  water  in  the  membrane,  nd  also 
depends  on  the  membrane  water  content.  According  to  Za- 
wodzinkski  et  al.  [23],  the  drag  coefficient  is  equal  to  unity  in 
the  range  of  water  content  from  A  =  5  to  14.  Because  the  drag 
coefficient  must  drop  to  zero  at  A  =  0,  below  a  water  content 
of  5,  a  linear  relation  between  the  drag  coefficient  and  wa¬ 
ter  content  is  used.  The  drag  coefficient  for  the  membrane  in 
equilibrium  with  liquid  water  increases  from  nd  =  1  to  2.5  as 
the  water  content,  A  increases  from  14  to  22  reported  by  the 
previous  study  of  Zawodzinkski  et  al.  [24]. 


n  d 


f  0.2A,  for  A  <  5 

1,  for  5  <  A  <  14 

l  0.1875A  -  1.625,  for  A  >  14 


(14) 


After  the  water  concentration  is  calculated  in  the  membrane 
with  Eqs.  (11)  and  (12),  the  water  content  of  the  membrane 
is  calculated  with  [4] : 


EWm  •  cw 
Pm  *  ^H20 


(15) 


where  pm  is  the  dry  membrane  density  and  EWm  is  the  equiv¬ 
alent  molecular  weight  of  the  membrane. 


3.3.  Boundary  conditions 

The  governing  equations,  boundary  conditions  and  normal 
vectors  at  the  inside  boundaries  are  summarized  in  Fig.  2,  the 
other  remaining  boundaries  are  either  insulation  or  symmetry. 
Due  to  symmetry  only  the  half  of  the  gas  channels  and  bipolar 
plate  shoulders  are  modeled. 


where  Nw  is  the  net  water  mass  flux  vector.  Positive  values  for 
Nw  mean  net  water  flux  from  anode  to  cathode  and  negative 
values  are  from  cathode  to  anode. 


3.3.1.  Catalyst  reactive  boundaries 

Water  can  enter  the  membrane  only  from  the  anode  and 
cathode  boundaries.  The  pressures  at  these  boundaries  are  set 
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to  the  solved  for  anode  and  cathode  pressures.  However,  the 
boundary  condition  for  the  water  concentration  offers  a  chal¬ 
lenge  due  to  the  hydrophobic  nature  of  the  membrane.  Over 
the  entire  range  of  normal  operating  conditions,  the  activity 
coefficient  of  water  in  the  membrane  is  greater  than  unity  if 
one  assumes  a  Raoult’s  law  relationship  between  water  ac¬ 
tivity  and  membrane  water  content.  In  theory  one  expects 
that  the  membrane  in  contact  with  saturated  vapor  and  liquid 
water  to  have  the  same  amount  of  water  content,  however 
experiments  conducted  on  several  polymer  polymer/solvent 
systems  reports  different  water  uptake  from  liquid  water  ver¬ 
sus  saturated  water  vapor.  The  phenomena  was  first  reported 
in  1903  by  Schroeder  and  thus  called  Schroeder’s  paradox. 

The  hydrophobic  nature  of  the  membrane  and  Schroeder’s 
paradox  are  treated  with  equivalent  boundary  conditions.  In 
this  approach,  a  fit  of  the  experimental  relationship  of  the 
water  content  of  the  membrane  to  the  water  vapor  activity,  a 
reported  by  Springer  et  al.  [2]  is  used  to  calculate  the  water 
content  at  the  membrane  boundaries  using: 


Table  2 


Butler- Volmer  kinetic  parameters 


Symbol 

Value 

Reference 

af,  anode;  anodic  transfer  coefficient 

0.5 

[13] 

af,  anode;  cathodic  transfer  coefficient 

0.5 

[13] 

af,  cathode;  anodic  transfer  coefficient 

1 

[13] 

cathode;  cathodic  transfer  coefficient 

1 

[13] 

ka,  anode  rate  constant  ((Am-2)  (m3  mol-1)0  5) 

533 

Estimated 

kc,  cathode  rate  constant  ((Am-2)  (m3  mol-1)) 

0.018 

Estimated 

yo0 ,  concentration  parameter 

1 

[8] 

Xh-,  ,  concentration  parameter 

0.5 

[8] 

with  the  following  expressions  [1] 


\  m2  /  \  XH+,  a 

,  _  -ref  /  £H2  \  /  fH+ 

0,a  0,a  |  ^ref  I  I  ref 

-H2  /  V  CH+ 


^0,c  — 


kH+ 


,C 


(22) 

(23) 


X  =  0.043  +  17.81a  —  39.85a2  +  36a3,  forO  <  a  <  1 

(16) 


which  is  used  with  Eq.  (15)  to  calculate  the  concentration  of 
water  at  the  anode  and  cathode  gas  channel  boundaries. 

As  the  water  mole  fraction  exceeds  saturation,  a  linear 
relation  is  assumed  between  the  water  content  and  water  ac¬ 
tivity  (Springer  et  al.  [2]): 


A  =  14  +  1.4(a  —  1),  for  1  <  a  <  3 


The  activity  in  the  vapor  phase  is: 

th2o  P 
a  = - 


P  sat 


(17) 


(18) 


where  psat  is  the  saturation  pressure  of  water  calculated  with 
[25]: 


where  i and  are  the  reference  exchange  current  densi¬ 
ties  for  anode  and  cathode  reactions  at  reference  conditions, 
and  y  is  the  reaction  orders  with  respect  to  hydrogen,  proton, 
and  oxygen  for  the  anode  and  cathode  reactions.  Assuming 
concentration  of  protons  at  the  anode  and  cathode  is  constant, 
Eqs.  (22)  and  (23)  can  then  be  simplified: 

*0,a  =  &a  •  Cu2m2  (24) 

io,c  =  kc  •  co2y°2  (25) 

where  ka  and  kc  are  constants  that  dependent  on  the  catalyst 
layer  properties.  Table  2  lists  the  known  parameter  values  of 
Eqs.  (20)-(25). 

For  the  cathode  reaction,  the  activation  overpotential  is 
defined  as: 


0  —  0s  0m  Voc 


(26) 


P  sat  —  exp 


f  73.648 


7258.2 

T 


7.3037  logT  +  4.1653  x  10~6T2) 

(19) 


Because  the  catalyst  layers  are  very  thin  compared  to  other 
elements  of  the  fuel  cell,  they  are  treated  as  reactive  bound¬ 
aries.  The  Butler- Volmer  kinetic  equation  is  used  to  obtain 
the  local  current  density  distribution  at  the  catalyst  surface 
[26]. 


—  ^0,a 


exP  1  aa  I  ~  exP  [  ~a 


h  —  ^0,c 


exp  I 


F 

~RT 


r]  —  exp  —a 


_  F 

c 

+ F 


0 


RT 


T] 


(20) 


(21) 


where  z'o,a  and  z'o,c  are  the  exchange  current  densities  of  the  an¬ 
ode  and  cathode,  respectively,  r)  the  activation  overpotential, 
and  a  is  the  transfer  coefficient.  The  exchange  current  densi¬ 
ties  are  scaled  from  the  reference  exchange  current  densities 


where  Voc  is  the  thermodynamic  open  circuit  potential  for  the 
overall  reaction: 


02  +  2H2  ->  2H20  (R-3) 

and  is  calculated  using  the  Nernst  law  [26], 

foe =1-229  -  0.9  X  10_3(7’  -  298)+2.3^  log(pl  po2) 

4 1  z 

(27) 


This  reduces  to  [27]: 


Voc  =  0.2329  +  0.0025  x  T  (28) 

For  the  anode  reaction,  the  activation  overpotential  is  the  dif¬ 
ference  between  the  electrode  and  the  membrane  potentials 
and  the  open  circuit  potential  of  the  anode  is  based  on  the 
standard  hydrogen  electrode: 

0  =  0s  0m 


(29) 
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The  concentrations  at  the  reactive  boundary  are  calculated 
from  the  mole  fractions  of  species  obtained  by  the  transport 
relations  shown  in  Fig.  2. 

For  the  catalyst  layer  between  the  anode  GDE  and  the 
membrane,  the  current  fluxes  are  given  with  Eqs.  (30)  and 

(31): 

n  •  (— crsV0s)  =  -/a  (30) 

n  •  (-<rmV0m)  =  4  (31) 

The  permeability  of  the  membrane  for  all  components  except 
water  is  assumed  to  be  zero  and  the  consumption  of  reactants 
and  formation  of  products  at  the  catalyst  reactive  boundary 
are  obtained  by  coupling  the  hydrogen  flux  to  the  local  current 
density  at  the  boundary. 


Table  3 


Physical  property  parameters 


Symbol 

Value 

Reference 

k,  permeability  of  the  GDE  (m2) 

1.76  x  1CT11 

[3] 

/ra,  viscosity  of  anode  gas  (Pas) 

1.378  x  10“5(r/298)102 

[32] 

Uc,  viscosity  of  anode  gas  (Pas) 

1.094  x  10“5(r/298)105 

[32] 

Ui,  viscosity  of  water  in  the 

3.56  x  10“4 

[1] 

membrane  (Pa  s) 

of ff,  conductivity  of  the  GDE 

570 

[31] 

(S  m-1) 

kp,  hydraulic  permeability  of 

1.8  x  10"18 

[3] 

membrane  (m2) 

pnf ,  dry  membrane  density 

1980 

[10] 

(kgm  3) 

EWm,  equivalent  molecular 

1.1 

[10] 

weight  (kg  mol  1 ) 


n  ■  I  -pwUl  (  D^j'Vxj 

V  ;=H2,H20  ^ 

Vi 

+{Xj  -  Wj ) — 

P 

The  mass  flux  of  water  at  the  anode  GDE-cataly  st-membrane 
boundary  is  equal  to  the  membrane  water  flux  at  the  boundary 
and  overall  mass  balance  of  water  at  anode  is  obtained  by 
coupling  through  the  momentum  balance  with  the  following 
equation: 


+  pwu2u 


l  Q 

= - M\  \- 

2  F  - 


(32) 


(-j^Mu2  -  Nw  j 

n  u  =  ± - L  (33) 

P 

The  current  flux  at  the  cathode  boundary  between  electrode 
and  membrane  are  given  with  Eqs.  (34)  and  (35),  respectively. 


n  •  (-<tsV0s)  =  4  (34) 

n  •  (  cym  V (pm)  =  4  (35) 


The  mass  flux  of  oxygen  at  the  cathode  catalyst  boundary  is 
obtained  by: 


n  ■  |  -pwo2  ^3  (  D°2 i^xi 

j= C>2,H20,N2 


V  p\  \  ic 

+(xj  -  Wj)— j  +  pwo2u  J  =  — Mo 


For  water: 


n  •  |  -pwu2o  ^  (  Dh2o j^Xj 

j=02,H20,N2 


(36) 


And  the  momentum  flux  at  the  boundary  is  given  by: 

o2  ^4^h2o  +  ^Vw) 

n  u  =  A - L  (38) 

P 

The  conservation  of  nitrogen  in  the  cathode  is  obtained  by 
solving  Eqs.  (36)-(38)  together,  yielding  zero  net  nitrogen 
flux  at  the  GDE-catalyst-membrane  boundary. 

3.3.2.  Bipolar  plate  shoulder-GDE  boundaries 

The  mass  and  momentum  transport  boundary  conditions 
between  the  bipolar  plate  shoulders  and  the  gas  distribution 
electrodes  are  all  insulation.  The  boundary  condition  for  the 
charge  balance  in  the  electrodes  is  set  to  zero  potential  at  the 
anode  and  set  equal  to  the  operating  cell  potential  Vceii  at  the 
cathode. 

3.3.3.  Channel-GDE  boundaries 

The  fuel  and  oxidant  streams  enter  and  leave  the  GDEs 
through  the  boundary  between  the  channels  and  GDEs.  The 
pressure  of  the  fuel  and  oxidant  streams  are  set  at  these  bound¬ 
aries  to  the  anode  and  cathode  operating  pressures  times  the 
mole  fraction  of  components.  Pure  hydrogen  is  used  as  a  fuel 
and  air  is  used  for  the  oxidant.  The  fuel  and  oxidant  streams 
are  fully  humidified  at  the  cell  operating  temperature  and 
pressures  for  both  anode  and  cathode  for  the  base  case. 

The  partial  pressures  of  water  at  the  anode  and  cathode 
boundaries  are  calculated  from  Eq.  (19),  and  are  divided  by 
the  operating  pressure  to  obtain  the  mole  fraction  of  water. 
The  remaining  components  are  hydrogen  at  the  anode,  and 
oxygen  and  nitrogen  with  a  stoichometric  ratio  of  0.21/0.79 
at  the  cathode. 

The  physical  property  parameters  are  shown  in  Table  3. 


4.  Solution  technique 


+  pwu2ou 


) 


4 

-  — 37ii2o  +  Nw 
Zr 

(37) 


A  finite  element  computational  fluid  dynamics  package, 
FEMLAB®  was  used  to  solve  the  non-linear  system  of 
equations  [19].  The  following  physics  application  modes  in 
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FEMLAB®  and  FEMLAB®  Chemical  Engineering  Module 
were  used  for  solving  the  dependent  variables: 

•  Darcy’s  law  for  pressure  in  the  porous  GDEs  ( p)\ 

•  Maxwell-Stefan  multi-component  diffusion  and  convec¬ 
tion:  for  weight  fraction  of  oxygen  and  water  at  the  cath¬ 
ode  gas  channels  and  GDE  ( wq2  and  u;h2o,c)  and  weight 
fraction  of  hydrogen  at  the  anode  gas  channels  and  GDE 

Oh2); 

•  conductive  media  DC  application  mode  for  charge  balance 
in  GDE  and  membranes  (0S  and  0m); 

•  general  PDE  form  is  used  to  solve  for  hydraulic  pressure 
and  water  concentration  in  the  membrane  (pw  and  cw). 

The  equations  and  boundary  conditions  were  outlined  in 
Fig.  2. 

FEMLAB®  uses  a  triangular  mesh  for  2D  geometries.  Ex¬ 
tensive  numerical  tests  were  performed  and  it  was  found  that 
1 800  mesh  elements  provided  satisfactory  spatial  resolution 
for  the  base  case  geometry  (Table  4)  and  the  solution  was 
found  to  be  independent  of  the  grid  size  with  further  refine¬ 
ment. 

A  stationary  non-linear  solver  was  used  together  with  Di¬ 
rect  (UMFPACK)  linear  system  solver.  The  relative  tolerance 
for  the  error  criteria  was  1  x  10~4  and  because  the  depen¬ 
dent  variables  vary  greatly  in  magnitude,  manual  scaling  of 
the  dependent  variables  was  used  to  improve  numerical  con¬ 
vergence.  The  manual  scaling  values  are  kept  constant  and 
were  selected  such  that  the  magnitude  of  the  scaled  degrees 
of  freedom  was  equal  to  one. 

The  model  was  created  and  tested  in  FEMLAB®  graph¬ 
ical  user  interface  and  then  saved  as  a  MATLAB®  M-File 
[28].  Performance  curves  were  obtained  by  varying  the  oper¬ 
ating  potential,  i.e.,  Vceii  within  a  loop.  Sensitivity  analysis  of 
cell  design  and  operating  conditions  on  the  fuel  cell  perfor¬ 
mance  was  conducted  by  putting  “for  loops”  in  MATLAB® 
for  changing  parameters  and  calling  the  FEMLAB®  non¬ 
linear  solver  with  updated  parameters.  The  model  geometry 
was  updated  and  re-meshed  only  when  a  geometric  dimension 
was  changed  thus  saving  computational  time.  The  2D  simula¬ 
tion  for  each  operating  cell  potential  converged  in  30-400  s, 
the  larger  being  for  the  high  current  density  and  limiting  re¬ 
actant  cases.  The  majority  of  the  runs  were  completed  under 
150  s  on  an  Intel  Pentium®  4  3.2  GHz  CPU  with  1  GB  of 
DDRam. 


5.  Results  and  discussion 

5.1.  Model  validation 

The  anode  and  cathode  reaction  constants,  ka  and  kc  were 
adjusted  to  fit  the  published  single  cell  experimental  results 
of  Ticianelli  et  al.  [29].  The  average  current  density  at  the 
catalyst  layer  was  calculated  with  Eq.  (39),  and  compared 


Fig.  3.  Comparison  of  the  predicted  and  measured  [29]  cell  polarization 
curves. 


Table  4 


Base  case  geometric  parameters 


Parameter 

Value 

Lc h,  half  of  the  gas  channel  length  due 
to  symmetry  (m) 

Lsh,  half  of  the  bipolar  shoulder  length 
due  to  symmetry  (m) 
ta  and  tc,  anode  and  cathode  gas  dis¬ 
tribution  electrode  thickness  (m) 
tm,  membrane  thickness  (m) 
s,  gas  porosity  of  the  anode  and  cath¬ 
ode  GDE 

0.5  x  10“3 

0.5  x  10“3 

2.6  x  10“4 

2.3  x  10“4 
0.6 

with  experimental  results  as  shown  in  Fig.  3. 

^  fl^sh  5~Lch 

4vg  —  T  T  /  i(y)dy 

^sh  T  Pch  Jo 

(39) 

The  fuel  cell  dimensions  and  operating  conditions  are  tabu¬ 
lated  in  Tables  4  and  5,  respectively.  The  estimated  reference 
exchange  current  density  for  the  anode  is  so  large  that  the  ac¬ 
tivation  overpotential  is  very  small  compared  to  the  cathode 
activation  overpotential,  which  is  in  accordance  with  experi¬ 
mental  findings  [29] . 

Table  5 

Operating  condition  parameters 

Symbol 

Value 

T,  temperature  (K) 

pa,  anode  side  pressure  (Pa) 

pc,  cathode  side  pressure  (Pa) 

Vq2  ,  cathode  feed  oxygen  mole  fraction 

C’  cathode  feed  water  mole  fraction 
,  cathode  feed  nitrogen  mole  fraction 
,  anode  feed  hydrogen  mole  fraction 
v^2o  a,  anode  feed  water  mole  fraction 

353 

303975 

506625 

0.1904 

0.0934 

i  —  x(A  —  4  n  „ 

02  H20,C 

0.8444 

1  ~xh2 
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Fig.  4.  Membrane  water  content  (A)  and  net  water  flux  vectors  at 
three  different  current  densities:  (a)  Vceii  =0.9  V,  i.e.,  /avg  =  0.05  A  cm-2, 
(b)  yceii  =  0.7  V,  i.e.,  /avg  =  0.52  A  cm-2,  and  (c)  Vceii  =  0.5  V,  i.e., 
/avg  =  1.13  A  cm-2  (left  side  is  the  anode  and  right  side  is  the  cathode). 


5.2.  Water  balance 

While  the  published  experimental  current  voltage  curves 
give  data  on  performance  of  the  cell,  membrane  water  content 
data  is  limited  [30].  The  calculated  membrane  water  content 
and  the  net  water  flux  are  shown  in  Fig.  4  for  three  differ¬ 
ent  operating  potentials  (current  densities).  At  low  current 
densities  the  water  diffusive  flux  is  from  anode  to  cathode 
as  the  water  content  is  higher  at  the  anode  side  in  Fig.  4a. 
At  low  current  densities  the  migration  of  water  due  to  proton 
flow  and  the  diffusive  flux  is  less  than  the  convective  flux, 
thus  the  net  water  flux  is  from  cathode  to  anode  shown  by 
the  net  water  flux  vectors.  As  the  current  density  increases 
the  amount  of  water  generated  at  the  cathode  side  increases 
and  thus  the  water  content  of  the  membrane  is  higher  at  the 
cathode  side.  Even  though  the  diffusive  water  flux  is  from 
cathode  to  anode,  the  migration  flux  overcomes  the  diffusive 
and  convective  fluxes  and  the  net  water  flux  is  from  anode 
to  cathode  (Fig.  4b).  The  convective  flux  of  water  is  from 
cathode  to  anode  due  to  the  pressure  gradient  between  the 
two  sides  and  does  not  change  with  increasing  current  den¬ 
sity.  The  effect  of  the  bipolar  shoulders  on  the  water  content 
of  the  membrane  can  clearly  be  seen  in  Fig.  4b  contours. 
The  water  content  gradient  is  higher  at  the  upper  section  of 
the  model  geometry  where  the  bipolar  plates  touch  the  elec¬ 
trodes.  As  the  cell  is  operated  at  higher  current  densities  the 
migrative  flux  exceeds  the  diffusive  and  convective  fluxes  sig¬ 
nificantly  and  dehydration  of  the  membrane  at  the  anode  side 
where  the  bipolar  plate  touches  the  electrodes  can  be  seen 
in  Fig.  4c. 

The  average  water  content  of  the  membrane  was  calculated 
by  integrating  the  water  content  in  the  membrane  and  dividing 


Fig.  5.  Effect  of  channel  geometry  on  the  cell  performance. 


it  by  the  total  membrane  area: 

^  /‘Z'sh  +  Z'ch  ptm 

A  avg  =  — - — — — —  /  /  Ux,  y)dx  •  d  v  (40) 

(^sh  T  Pch)  ’  hn  JO  JO 

At  low  current  densities  Fig.  4a,  the  average  water  content 
of  the  membrane  was  12.57  mol  H2O  (mol  SO3-)-1  and  as 
the  current  density  increases  the  average  water  content  in  the 
membrane  reaches  13.83  in  Fig.  4b  and  drops  to  13.06  at  high 
current  density  case  Fig.  4c. 

5.3.  Channel  size  effects 

Sensitivity  analysis  of  the  channel  dimensions,  in  the  range 
of  1-2.5  mm,  on  the  fuel  cell  performance  is  shown  in  Fig.  5. 
In  all  of  these  cases,  the  ratio  of  gas  distribution  channel  to 
bipolar  plate  shoulder  was  kept  constant  at  1 ,  i.e.,  equal  chan¬ 
nel  and  bipolar  shoulder  size.  The  effect  of  channel  dimen¬ 
sions  on  the  fuel  cell  power  density  becomes  more  important 
at  high  current  density  applications.  For  example,  a  1  mm 
channel  and  bipolar  shoulder  has  14%  more  power  than  the 
2.5  mm  channel  and  shoulder  design  while  operated  at  0.6  V; 
however,  power  density  difference  is  11%  while  the  cell  is 
operated  at  0.7  V.  Performance  curves  for  the  2.5  mm  channel 
and  shoulder  case  clearly  shows  mass  transfer  limitations  at 
current  density  above  0.8  A  cm-2.  This  is  shown  by  a  drop 
in  the  performance.  The  model  predicts  that  a  current  den¬ 
sity  of  1  A  cm-2  cannot  be  obtained  with  2.5  mm  design  due 
to  mass  transfer  limitations.  The  1  mm  channel  and  shoulder 
case  shows  a  constant  slope  for  moderate  to  high  current  den¬ 
sity  operations  up  to  1 .4  A  cm-2  and  it  can  be  concluded  that 
mass  transfer  limitations  to  catalyst  sites  under  the  bipolar 
plate  shoulder  can  be  prevented  by  using  smaller  channels. 

5.3.1.  Current  density  profiles 

Fig.  6  shows  the  current  density  profiles  at  the  cathode 
catalyst  layer  for  the  1  and  2  mm  channel  and  shoulder  cases 
at  three  different  operating  cell  potentials  0.9,  0.7,  and  0.5  V. 


104 


G.H.  Guvelioglu,  H.G.  St  enger  /  Journal  of  Power  Sources  147  (2005)  95-106 


Fig.  6.  Current  density  profiles  at  the  catalyst  layer  for:  1  mm  channel  and 
1  mm  shoulder,  2  mm  channel  and  2  mm  shoulder. 


Fig.  7.  Hydrogen  mole  fraction  profiles  at  the  anode  catalyst  layer  for:  1  mm 
channel  and  1  mm  shoulder,  2  mm  channel  and  2  mm  shoulder. 


This  figure  shows  low,  moderate,  and  high  current  density 
operations,  respectively.  The  channel  and  bipolar  plate  was 
simplified  with  symmetry  at  the  center  sections  of  each  region 
so  only  0.5  mm  of  the  channel  and  0.5  mm  of  the  shoulder 
is  shown  in  the  plot  for  1  mm  case  and  1  mm  of  channel  and 
shoulder  is  shown  of  2  mm  case.  Length  =  0  is  the  center  of  the 
gas  channel  and  at  length  =  0.5  mm  the  bipolar  plate  shoulder 
starts  touching  the  electrode  for  the  1  mm  case.  At  low  cur¬ 
rent  density  operation,  i.e.,  Well  =  0.9  V,  the  current  density 
profile  is  uniform  for  both  cases.  At  moderate  current  density 
operation,  i.e.,  Vceii  =  0.7  V  the  current  density  increases  as 
it  gets  closer  to  the  location  where  the  bipolar  plate  touches 
the  electrode.  This  increase  is  due  to  the  low  conductivity 
of  the  electrode  as  well  as  increasing  water  concentration  of 
the  membrane  and  reduced  local  membrane  resistance.  The 
electrode  conductivity  of  570  S  m-1  was  estimated  from  the 
commercial  electrode  properties  and  in  agreement  with  the 
estimate  of  Nguyen  et  al.  [31].  As  we  pass  the  location  where 
the  bipolar  plate  shoulder  touches  the  electrodes,  a  uniform 
current  density  profile  is  observed  suggesting  that  the  concen¬ 
tration  change  of  reactants  does  not  effect  the  current  density 
profile  due  to  also  decreasing  membrane  resistance.  When 
the  current  density  profiles  at  the  high  current  density  opera¬ 
tions  was  studied,  i.e.,  Vceii  =  0.5  V,  the  current  density  profile 
drop  was  significant  due  to  a  concentration  drop  of  reactants. 

5.3.2.  Concentration  profiles 

Figs.  7  and  8  show  the  hydrogen  mole  fraction  at  the  anode 
catalyst  layer  and  oxygen  mole  fraction  at  the  cathode  catalyst 
layers,  respectively. 

At  low  current  densities  the  hydrogen  mole  fraction  pro¬ 
file  at  the  anode  catalyst  layer  shows  a  decrease  as  it  goes 
under  the  bipolar  plate  shoulder.  However,  at  moderate  and 
high  current  densities  the  hydrogen  mole  fraction  increases 
even  though  hydrogen  is  consumed  at  the  catalyst.  This  in¬ 
crease  is  due  to  the  complex  water  balance  in  the  PEM  fuel 
cell  as  shown  earlier  in  Fig.  4.  At  low  current  densities,  the 


net  water  flux  is  from  cathode  to  anode  and  the  hydrogen  is 
consumed  at  the  catalyst  layers.  This  results  in  a  drop  of  hy¬ 
drogen  mole  fraction  in  Fig.  7;  however,  at  moderate  to  high 
current  densities  the  net  water  flux  is  from  anode  to  cath¬ 
ode  in  the  membrane.  Even  at  moderate  current  densities,  the 
magnitude  of  the  net  water  flux  is  more  than  the  hydrogen 
consumption  rate  thus  yielding  an  increase  in  mole  fraction 
of  hydrogen  (Fig.  7  for  Vcen  =  0.7  V).  The  difference  of  hy¬ 
drogen  mole  fractions  at  length  =  0,  i.e.,  the  catalyst  layer 
below  the  center  of  the  gas  distribution  channel,  was  due  to 
the  difference  in  current  densities  as  shown  in  Fig.  6. 

The  effect  of  channel  and  shoulder  dimensions  on  the  oxy¬ 
gen  mole  fraction  at  the  cathode  catalyst  layer  was  more  sig¬ 
nificant  and  is  shown  in  Fig.  8.  The  complex  water  trans¬ 
port  phenomena  in  the  membrane  as  well  as  consumption  of 
oxygen  magnifies  the  mole  fraction  drop  of  oxygen  at  high 
current  density  cases.  For  the  2  mm  channel  and  shoulder 
case  the  oxygen  mole  fraction  drops  to  0.03.  This  is  the  main 


Fig.  8.  Oxygen  mole  fraction  profile  at  the  cathode  catalyst  layer  for:  1  mm 
channel  and  1  mm  shoulder,  2  mm  channel  and  2  mm  shoulder. 
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Fig.  9.  Average  membrane  water  content  for  cell. 

reason  for  the  drop  in  current  density  under  the  shoulders  in 
Fig.  6  for  Vceii  =  0.5.  Also,  it  is  worth  mentioning  that  the 
2  mm  channel  and  shoulder  might  not  be  a  proper  design  for 
high  current  density  applications.  At  the  exit  of  the  fuel  cell, 
the  oxygen  mole  fraction  will  be  far  less  than  that  of  this 
case  and  may  lead  to  extinction  of  reactants  under  the  bipolar 
plates. 

5.3.3.  Water  content 

Channel  and  bipolar  plate  shoulder  dimensions  do  not  only 
effect  the  mass  transfer  limitations  but  also  effect  the  water 
content  of  the  membrane.  Fig.  9  shows  the  effect  of  the  chan¬ 
nel  geometry  on  the  average  water  content  of  the  membrane 
at  different  current  densities,  calculated  using  Eq.  (40). 

At  current  densities  below  0.2  A  cm-2  the  average  mem¬ 
brane  water  content  is  almost  the  same  for  all  geometries.  As 
the  current  density  increases,  the  membrane  water  content 
also  increases.  This  increase  is  due  to  the  increasing  concen¬ 
tration  of  water  at  the  cathode  catalyst  membrane  boundary. 
As  the  current  density  increases  above  0.2  A  cm-2  the  mem¬ 
brane  water  content  also  increases  but  the  rate  of  increase 
is  smaller  compared  to  current  densities  below  0.2  A  cm-2. 
This  is  due  to  an  increase  in  the  migrative  flux  from  anode  to 
cathode.  The  effect  of  channel  and  shoulder  dimensions  on 
the  average  membrane  content  can  be  seen  in  Fig.  9.  Larger 
channel  and  shoulders  favor  higher  membrane  hydration  at 
moderate  current  densities,  i.e.,  0.2-0. 8  A  cm-2  as  the  large 
shoulder  area  limits  the  water  transport  from  inside  the  elec¬ 
trodes  under  the  bipolar  plate  to  channels.  The  membrane 
water  content  reaches  its  maximum  value  for  all  channel  and 
shoulder  dimensions  between  0.4  and  0.5  A  cm-2  where  the 
net  flux  of  water  becomes  zero.  As  the  current  density  in¬ 
creases  above  0.8  A  cm-2,  smaller  channel  and  shoulder  di¬ 
mensions  favor  higher  membrane  hydration  levels.  As  the 
current  density  increases,  the  anode  side  of  the  membrane 
dries  as  can  be  seen  in  Fig.  4c,  thus  increasing  the  resistance 
of  the  membrane. 


Fig.  10.  Effect  of  relative  humidity  of  the  feed  on  the  performance  of  the 
cell. 

5.4.  Gas  distribution  electrode  porosity 

The  effect  of  porosity  of  the  gas  distribution  electrodes  was 
also  studied  and  it  was  found  that  a  reduction  of  porosity  from 
0.6  to  0.4  resulted  in  2.4%  decrease  in  the  average  current 
density  for  the  1  mm  channel  and  shoulder  case  operated  at 
0.6  V.  The  porosity  effect  was  found  to  be  more  significant  at 
higher  current  density  operations  and  a  6%  drop  was  observed 
for  a  fuel  cell  operating  at  0.4  V.  The  porosity  of  electrodes  is 
critical  for  cells  with  bigger  catalyst  regions  under  the  bipolar 
plate  shoulders  as  the  supply  of  the  reactants  can  be  rate 
limiting.  For  the  2  mm  channel  and  2  mm  shoulder  case  the 
reduction  of  electrode  porosity  from  0.6  to  0.4  resulted  in  a 
6.6%  drop  in  current  density  for  a  fuel  cell  operated  at  0.6  V. 
Due  to  insufficient  transport  of  reactants  to  the  catalyst  sites, 
the  2  mm  channel  and  shoulder  case  could  not  be  operated  at 
higher  current  densities  than  0.85  A  cm-2. 

The  water  content  of  the  fuel  and  oxidant  streams,  of  the 
anode  and  cathode  gases  effect  the  fuel  cell  performance  sig¬ 
nificantly  as  shown  in  Fig.  10.  The  effect  of  anode  gas  inlet 
relative  humidity  was  found  to  be  one  of  the  most  critical 
operating  conditions.  A  50%  relative  humidity  of  the  anode 
gases  was  not  enough  to  keep  the  membrane  well  hydrated 
thus  resistance  of  the  membrane  and  ohmic  loses  at  the  mem¬ 
brane  is  higher.  The  effect  of  cathode  relative  humidity  was 
not  as  critical  because  at  high  current  densities,  water  gen¬ 
erated  at  the  cathode,  and  the  flow  of  water  is  from  anode 
to  cathode.  At  low  current  densities,  less  than  0.4  A  cm-2, 
the  net  water  flux  was  from  cathode  to  anode  thus  the  rela¬ 
tive  humidity  effect  of  the  cathode  was  also  low  at  this  re¬ 
gion.  The  50%  relative  humidity  cathode  and  saturated  anode 
stream  performance  curve  drops  slightly  below  the  saturated 
anode  and  cathode  stream  performance  curve  between  0.2 
and  0.6  A  cm-2  regions.  As  the  cell  was  operated  at  higher 
current  densities  the  water  supplied  from  the  anode  side  as 
well  as  the  water  generated  at  the  cathode  was  sufficient  to 
keep  the  membrane  hydrated  at  high  levels. 
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6.  Conclusions 

A  computational  fluid  dynamics  model  of  a  polymer  elec¬ 
trolyte  membrane  fuel  cell  has  been  presented.  The  mass 
transport,  momentum  transport,  and  electrochemical  pro¬ 
cesses  occurring  in  the  membrane  electrolyte  and  catalyst 
layers  have  been  modeled.  The  effects  of  channel  and  bipolar 
plate  shoulder  size,  porosity  of  the  electrodes,  and  fuel  and 
oxidant  humidification  on  the  cell  performance  have  been 
studied  for  a  wide  range  of  current  densities. 

It  is  found  that  smaller  sized  channels  and  bipolar  plate 
shoulders  are  required  to  obtain  higher  current  densities, 
although  larger  channels  are  satisfactory  at  moderate  cur¬ 
rent  densities.  If  the  application  requires  bigger  channel  and 
shoulder  sizes,  increasing  the  porosity  of  the  gas  distribution 
electrode  helps  the  mass  transport. 

The  effect  of  the  relative  humidity  of  the  anode  gas  stream 
was  found  to  be  the  most  critical  condition  affecting  the  per¬ 
formance  of  the  fuel  cell.  Also,  the  fuel  cell  design,  the  ge¬ 
ometric  dimensions  of  the  channels  and  bipolar  plate  shoul¬ 
ders,  the  thickness  of  electrodes,  and  membrane,  and  the 
porosity  of  the  electrodes  and  conductivity  of  the  electrodes 
needs  to  be  selected  carefully  for  different  applications. 

The  development  of  this  FEMLAB®  model  enables  us 
to  study  operating  conditions  and  design  parameters  from 
a  CFD  model  and  optimization  algorithms  available  in 
MATLAB®.  The  selected  2D  domain  enabled  us  to  study  the 
effects  of  channel  and  shoulder  sizes,  which  along  the  channel 
2D  models  lack.  However,  the  concentration  changes  along 
the  channel  and  their  effect  on  the  performance  of  the  entire 
fuel  cell  cannot  be  captured  with  this  model.  Future  work  will 
focus  on  developing  a  procedure  to  evaluate  the  entire  fuel 
cell  performance  with  multiple  evaluations  of  the  2D  model 
along  the  channel.  The  operating  conditions  and  selection  of 
components  for  different  applications  can  then  be  optimized. 
The  robustness  and  typical  convergence  speed  of  under  150  s 
with  a  desktop  PC  is  a  crucial  benefit  of  this  model  and  will 
enable  CFD  optimization  to  be  a  practical  application. 
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